#include "config.h"

#include <clx/libECM.h>
#include <gsl/gsl_fit.h>
#include <gsl/gsl_vector.h>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_multifit_nlin.h>
#include <gsl/gsl_blas.h>

#include <math.h>

#include "lstatfit.h"

#define ERR(fmt, ...)	ECMERR("ERR: lstatarr: " fmt, ## __VA_ARGS__)

static char	*inputfile;

static ecm_uint32	n_arr_inputs;

static ecm_double	*iarrival_times;
static ecm_double	*probs;

static void
usage(void)
{
	ecm_fprintf(ecm_stderr,
"Usage: lstatarr [ <options> ] <arrival file>\n"
" Options:\n"
"   -h: help\n"
"   -o: output file\n"
	);
}

#define ALLOC_CHUNK	256

static ECMBOOL
read_arr_input(void)
{
	FILE	*fp;
	char	buf[256];
	ecm_uint32	n_alloced = ALLOC_CHUNK;

	fp = ecm_fopen(inputfile, "r");
	if (fp == NULL) {
		ERR("%s: cannot open file", inputfile);
		return ECM_FALSE;
	}

	iarrival_times = (ecm_double *)ecm_malloc(sizeof(ecm_double) * n_alloced);
	probs = (ecm_double *)ecm_malloc(sizeof(ecm_double) * n_alloced);

	while (ecm_fgets(buf, 256, fp)) {
		ecm_uint32	iarrival_time;
		ecm_double	prob;

		if (*buf == '#')
			continue;
		if (ecm_sscanf(buf, "%u %lf", &iarrival_time, &prob) != 2) {
			ERR("%s: invalid format", buf);
			continue;
		}

		if (n_arr_inputs >= n_alloced) {
			probs = (ecm_double *)ecm_realloc(probs, sizeof(ecm_double) * (n_alloced + ALLOC_CHUNK));
			iarrival_times = (ecm_double *)ecm_realloc(iarrival_times, sizeof(ecm_double) * (n_alloced + ALLOC_CHUNK));
			n_alloced += ALLOC_CHUNK;
		}

		iarrival_times[n_arr_inputs] = iarrival_time;
		probs[n_arr_inputs] = prob;
		n_arr_inputs++;
	}

	return ECM_TRUE;
}

static void
parse_args(int argc, char *argv[])
{
	int	c;

	while ((c = ecm_getopt(argc, argv, "h")) != -1) {
		switch (c) {
		case 'h':
			usage();
			ecm_exit(0);
		default:
			usage();
			ecm_exit(1);
		}
	}

	if (ecm_optind == argc) {
		usage();
		ecm_exit(1);
	}

	inputfile = ecm_strdup(argv[ecm_optind]);
}

int
main(int argc, char *argv[])
{
	msl_fit_t	*fit;

	ECM_initialize();

	ecm_set_output_file(ecm_stderr);
	ecm_set_verbose_mode(ECM_TRUE);

	parse_args(argc, argv);

	if (!read_arr_input())
		return 2;

	fit = lstat_fitting("arrival", MSL_FIT_EXPONENTIAL, n_arr_inputs, iarrival_times, probs);
	if (fit != NULL)
		lstat_report_fitting(fit, 0);
	else {
		ERR("cannot fitting arrival");
	}

	return 0;
}
